#include <stdio.h>
#include "septools.h"
    
int main(int argc, char **argv){
  data input;
  int file_counter=0; 
  double **inertia=sep_matrix(3,3);

  if ( argc != 2 ){
    fprintf(stderr, "Arguments: <type> \n"); 
    exit(EXIT_FAILURE);
  }

  char type = argv[1][0];
  get_num_type(&input, type, "00001_mol.conf");
  get_time(&input, "00001_mol.conf");

  allocate_memory(&input);

  FILE *fout = fopen("pmi.dat", "w");
  if ( fout == NULL )
    sep_error("Error openeing output file\n");

  double dt = input.tnow;
  while (1) {

    printf("\r");
    printf("Processing sample %d  ", file_counter);
    fflush(stdout);

    char file[256];
    sprintf(file, "%05d_mol.conf", file_counter);

    if ( file_exists(file) == 1 ) {
      read_entry(&input, file);

      double sum_eig[3] = {0.0}, eig[3];
      for ( unsigned i=0; i<input.npart_type; i++ ){
	inertia[0][0] = input.inertia[i][0];
	inertia[1][1] = input.inertia[i][1];
	inertia[2][2] = input.inertia[i][2];
	inertia[0][1] = input.inertia[i][3];
	inertia[0][2] = input.inertia[i][4];
	inertia[1][2] = input.inertia[i][5];
	inertia[1][0] = inertia[0][1];
	inertia[2][0] = inertia[0][2];
	inertia[2][1] = inertia[1][2];
      
	sep_eig_real_symmetric(eig, inertia);
	
	for ( int k=0; k<3; k++ )
	  sum_eig[k] += eig[k];
      }
      
      inertia[0][0] = input.inertia[0][0];
      inertia[1][1] = input.inertia[0][1];
      inertia[2][2] = input.inertia[0][2];
      inertia[0][1] = input.inertia[0][3];
      inertia[0][2] = input.inertia[0][4];
      inertia[1][2] = input.inertia[0][5];
      inertia[1][0] = inertia[0][1];
      inertia[2][0] = inertia[0][2];
      inertia[2][1] = inertia[1][2];
	
      sep_eig_real_symmetric(eig, inertia);

      fprintf(fout, "%f %f %f %f %f %f %f\n", dt*file_counter, 
	      sum_eig[0]/input.npart_type, sum_eig[1]/input.npart_type, 
	      sum_eig[2]/input.npart_type, eig[0], eig[1], eig[2]);

      file_counter++;
    }
    else {
      break;
    }

  }
  printf("\n");

  fclose(fout);
  free_memory(&input);
  sep_free_matrix(inertia, 3);

  return 0;
}

